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urpose of this cciputer program is to conpute the dynamic preesura 
distrisution and related quantities of interest due to the diffraction of sin- 
usoidal head waves. The method of computation is based on slender-body theory. 
The theoretical analysis is based on the assumption that the ship is slender. 
In addition, it is assumed that the incident waves are of small amplitude and 
their wavelength is short relative to the ship length. In the next section 
> a brief summary of the thecretical analysis in order to facili- 
. understanding of the computer program. Letails of the theoretical 
wsls may be found in Reck (1979). In the follcwing sections, details of 


conerical technique and the computer program wil] ve discussed. 


‘Tne coordinate syster is shown in Figure 1. The origin is at the bow 


with tne x-axis pointing aft. The z-axis is vertical upward and the x- lane 
y P 


waves 


Figure 1. Coordinate System 


is coincident with the caln-water plane. In this coordinate system, the 


incident wave potential is given by 


iw,t 
$, (x-yez)e © 


where ; 
Vz ~"1VX 
e e 


= acceleration gravity 
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= wave amplitude 
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= wave number 


= wi/g 
U = forward speed 


Wy = absolute wave frequency 
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The diffraction potential is written as 
iwet 
bp (x,y, ze He (2) 


To find ?p we must solve both the near-field and far-field problems. In the 
far-field, the solution is represented by a line distribution of pulsating 


sources, the strength of which is given by 
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The source strength, o(x), is found by solving the Volterra integral equaticr 
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kg (x) = near field source strength. 


Tne near-field source strength is found by solving the near field problem and 
May be assumed known when solving equation (3). 


In the near-field, the diffraction potential is written as 


a: pay ze (4) 
Beth che First and second order near-field solutions must be determined. The 
fire order sclution is simply the negative of the incident wave. The second- 
cide: solutio: ust satisfy the Nelmholtz eguetion subject to boundary condi- 
tiess on the free surface and the body surface, At infinity, the near-field 


eclucicon must match the far-field solution. 
As discussed in Beck (1979), the two-term near-field solution can be 
written as 


CR, 722) U7 a + atx] eB, (008 (¥.2) + ) By 2)Onty 2) 
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(5) 
9 coefficient A(x) is determined Ly matching witr the far-field solution and 
we find 


o (x) 


A(x) 3 = Tao tx) (6) 


The terms inside the square brackets represent Ursell's (1968a) solution 
to the Helmholtz equation subject to the free surface and body boundary condi- 
tions. $(y,z) is a source-like term and 0yn are wave-free potentials. The coef- 
ficients Bo, By are determined by satisfying the body boundary condition, which, 
in this case, is obtained by setting the normal derivatives of the terms in 
square brackets equal to zero on the body surface.? 

The determination of the coefficients Bo, By for arbitrary body shapes 
is often very difficult. Troesch (1976) has avoided this problem by using ar 


integral-equation technique. Troesch's technique and computer program were 


1. 7 1s type of solution is often called a multi-pole expansion. 
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actually developed for the solution of the oblique seas case, but they can be 
modified for use in the present problem. We first write Ursell's solution in 


the form 


+ 
Yly,Z2;x) =e de o(y,2;x) 
(7) 
The potential $(y,z,x) must satisfy the Helmholtz equation 
32 2 
Tua =< 
(8) 


subject to boundary conditions on the free surface and the body. At infinity, 
the behavior of ¥ must match the far-field solution. The free-surface boundary 


condition is 
—- vo =0 on 2=0 


On the body, the boundary condition is 


a Vz = 
oN [e + vtz,y20)| =0 


oN =- oN on h(y,z;x) = 0 (9) 


where 


h(y,z,x) = equation of body surface in y-x plane 
N = two-dimensional unit normal to body surface in the y-z plane 
= (No,N3) 


The positive sense of N is into the body. 
The problem for the potential % is now equivalent to the boundary value 
problem solved by Troesch. He writes the solution as a distribution over the 


body surface of two-dimensional sources. Thus, we have 


%(y,2:%x) =, | dl ¥(n,0)G(y,2:n,6) 


C(x) (10) 
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where the line integral is taxen along the bod; contour. y(n,¢) is the two-~ 


dimensional source streng7n, «nich is determine by satisfying the body boundary 


condition. G(y,23N, %) is <he Green function, wiich satisfies tne Helmholtz 


a 


equation and the free surZ2ce boundary condition (for details see Ursell 
(1968b)). By matching the limiting values of bo-vh the multipole solution and 


the Green function soluticn for large values of y, it can be shown that. 


, 
: 


Bo ix) = 2 | dl y(nite™ 
C(x) 


(12) 


Equation (11) allows the cetermination of the nenr-field source strength, Bo(x), 


by using the integral-equation solution technique without developing the multi- 
pole expansion. 

The pressure acting cn the body is found using the linearized Bernoulli 
equation, which may be writ<en as 

- iWet 

p= . 
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Substituting the expressicrs for $7 and ¢p into equation (12) and retainins 
only the lowest-order ter-s (see Beck (1979)), we find the following expression 


for the nondimensional pressure amplitude: 
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Likewise, we find the linearized wave amplitude in the near field is given 


by the expression 


S(x,y) — . o*(x) ; -ivx 
re a i TB (x) v(y,O;x)e 


(14) 


The exciting forces and moments are found by integrating the pressure over 
the body surface. As shown in Beck (1979), the vertical, sectional exciting 


force can be written as 


0 (x) ~ivx 


£(x)e% = E By (x) faery 200] e 
C(x) 


(15) 
The total heave exciting force and pitch moment about midship are giver. 


by the following integrals along the ship length: 


L “ 
F, = | dx f(x) ¢ =F 

? (16) 
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where 

F3 = total heave exciting force 


Fs = total pitch exciting moment 
L = ship length 
The wave induced bending moment is found by twice integrating the vertical 
exciting force up to the desired station. Setting the shear at the bow equal 


to zero and integrating by parts once, we arrive at 


x 
a f a 
BM(x) = -] a&(E-x) £,(6)e 
0 
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In the computer progran, the value of x, the station at which the bending 


nonent is evaluated, is specified by the user. 
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For any station whict. is a luplicate of the previous staticn (4.6. paralle. ~t1- 
body), tne results ef jas previous station are copied rathoe Uban resolving une 


near-fieia promen ‘WW. near-Tield proble~ is solved usii, the sethod of 


Troesch (1979). Because Ircesch's method has been accumented elsewhere (Tre=sen 


ig é mya s 4. oe } Ried a 
(19772). Froeech (2097° 4%. it “iil not be discus: 


ed here. His vompuber toccvan 
has mee. "ciiried to run in nead seas by eliminating the isaginary pur of ines 
Breen function and setting k=v. Many of Troesch's subroutines have been usec 
directiy and othersneeded cenly slight modification. I+ should be roted thet 
Troan ws acoee en won .eitten in double precision, so the’. cere cwitehtuge of 
mole: %. Wee. Lhe varintles in the present program ana the veriables in tre 
Trou.2) suivoutlines is necessary. This switching occurs in SUEROUY La SH lUr 


ana SUSROUTIDS TWDATA. 


For accuracy, the Volterra integral equation is sulved at 2cre sballens 
along the chin lengih chan just the input stations. The input ovations ers 


located at XAXIS(I). ‘Tne stations for the integral equation ani subseaue:- 
caleulaticis are locaisd at XI(I). The increased number of stations is developed 
in SUBROUTINE INSENT. nis subroutine inserts more ste lion. wo thal theis 


spacing is avprosimatcl. a cosine distribution alony the lengin. Phe ing: 


stations are never closer than t EPSIL to the input station. At present EPSIL 


is set eaual to .5% of the ship length. Figure 2 shows the two eis 


along the length. In ine figure, there are 7 input st.ticns ani 13 XI-stations. 
The array TSNUN(I) is used to store the number of the XI-station at each inzut 


station. For example, in Figure 2 ISNUN(4)=7. The values ef the various quantities 


input 
station 
number = 


Xf station { 11 12 13 
number 


Figure 2. Numbering of the XAXIS(I) and XI(TI) 
axis systems 


which are computed in the two-dimensional problem at the input station are 
determined at the XI-staztior by cubic-spline interpolation in SUBROUTIKE 
INTRPL. 

The Volterra integral equation (equation (3)) is solved at each XI station 
by a marching process starving at the bow. o(x) is assumed to vary linearly 
over each segment. The value of o(x) between the Ej] and §5 station is 
then give: by 


(E-£-) (544175) 
a noe es 7 sac 
5 j j 2 (19) 
where 
6: = value of o(x) at station number j 


0,5 = value of 7(x) at station number jt] 
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To develop the marching process we first rewrite the integral equation as 
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(20) 
where 30 541 equals the value of 39 at station number }+1. The integral in 
equation (20) can be evaluated analytically by substituting equation (19) as 
follows 


Xk+1 . : 
ae 2 | ve (ES) 4 * (Sq 4178) OK 
Ay VE = 
Xk jt+l ° 
Sk+]1 
= | g as 2m Ss. ’ 
Ay IG (Eres 4y re ote e541 5) 
(21) 
ok H : 
as ay IGS e541 Ob? = se hee as 
42 2, Jay 
where G(a,B,y) = (2a- a ae (22) 
Multiplying both sides cF ectation (20) by 6B Gy and substituting equation § 


(22) we arrive at the following expression for 96. 
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(23) 
In equation (23) the funccicn G(3,8,Y) is given by equation (22), 


To s.art the marching trocess, it is assumed tnat the nezr-field and thve:- 


dimensic:::2 source streng7.s sre zero at the bow (i.e. 0,=Bo, = 0). By 
expandir.; tne integral equzticn for o(x) around x=0 and taking the linit as x0, 
it can be shown that 


o'(x) = -1 ge B'(x) at x=0 
Ww ° 


0 (15) 
where the prime denotes diiferentiation with respect to x. The linear approxi- 
mation used for o(x) in the runerical scheme leads to exactly this same slope 
at x=0. 

At the stern, there exe *wo possible cases. For cruiser sterns, a(x) 
and Bo(x) are set equal tc zero. For transom sterns, the progra:: uses the 
values of Bo(L) and o(L) <: zonputed for the transom stern secticn. Thus, 
o(x) and Bo(x) are not ea21 to zero at the stern section. The validity of 
this result should be fircsner investigated, but it seems to give reasonable 
answers. 


T2 pressure and neaz-lield wave amplitude are computed by equations (13) 


and (1':.. In computing t>ss¢ quantities, the quotient o(x)/mBo(x) car not be 


computed fer points at which Bo(x)=0 (i.e. 2t the bow and for cruiser sterns). 
The proper limits for the quotient are found by expanding the integral equation 


for o(x) around x=0 or x=L and taking the proper limits. At the bow, we fird 


x70 oO 
At the stern, the linit need only be taken for cruiser type sterns. in this 


ease, the result is 


The exciting forces and midship bending moment are found from the evalua- 
tion of equations (15), (16), (17) and (18). These equations all involved 


integrals of the forr 
x 


I= | dé £(E)e2¥8 


0 
where £(&) is a function which varies smoothly along the length. To evaluaze 
these integrals a simple Filon-Trapezoidal rule is used (see Tuck (1967). This 


integration is carried out in SUBROUTINE TRAP. 


REQUIRED TRPUT 


The required input is subdivided into two separate parts. The first rari 
contains all the control information such as number of stations, wave frequencies, 
ship speed, etc. The second part is the offsets which describe the ship hull. 
The control information is read in on device number 5. The ship offsets are 


read in on device 7. 


Control Information (Read in on device 5) 


Card 1: FORMAT (715) 
NSTA - number of ship stations to be read in 


ND ~ nunber of divisions to be used in developing cosine station 
spacing used in setting up the XI(I) axis system. 


NVEL - number of velocities at which calculations are to be made 
(maxinun cf 8) 


NFREQ - number of wave frequencies at which calculations are to be made 
(maxinum of 16) 


number of stations at which pressure information is desired 
(maximun of NSTA) 


NPRES 
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MvAve = control constant for wave arplitrte caleuletic:. 
= One dave anplitude calculaziur 
=1 compute wave amplitude along ship sicz 


NRID  - the number of the station at which the midship bending monert is 
computed. 
Card 2: FORMAT (4F10.4:} 
RHO - water dunsi*.y (slugs/ft) 
GPAY - accelerstion of gravity (ft/sec) 


- ship length (7+) 


- incident wave splitude (ft) 


Care. 7% One (OF10.4) 
Frovn( 7) - Provue nunbers «4 which calculatien are to ae vade There snovwi 
be iVEL values of FRCUB(T). 
Car eer FCrUAT (GF1LO.E) 
WI) etd} Wavelength-to-ship-length ratios av which calculations are io 24 
nade. NESEQ values of WIXL(T) should be read in. 
Card 5: + GitiAT (1015) 
nUE(T) - The nunbers of the stations at which the pressure distributior 


is to be completed. There should be NPRES values of NUM(I). 
NOTE if NPRES=!STA or NPRES=0 data card 5 must be omitted. 


Ship Cffsets (Read in on device 7) 


There is one comp_ete set of data cards for each ship station. At leas- 
one offset (up to a maxinum of 25) must be given for each station. Only one 
offset should be given at the bow; it may be either the point (0,0) or, for e 
plum bow, the point (0,T) where T is the fore-foot draft. The ship may have 
a cruiser (one offset. point) stern or a transom (many offset points) stern. 

Figure 3? is a picture of the means by which the offsets are read in for a 


given station. For accuracy in the calculations, more points should be entered 


near the waterline and in areas where the shape changes rapidly (i.e. the tven 
of the bilse). The po‘n:zs are always entered starting at the negative waterline 
and reading counterclockwise. Only half the section should be entered. There 
must be data points at (-B/2,0) and (0,T). 


For each station, the following data cards are necessary. 


Card 1: FORMAT (15,F10.4) 


NT2(ISTA) = number of offset points on the hal? section for the station 
number ISTA. 


ave 
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YAXIS(5,1), 
ZAXIS (5,1) 


YAXIS(5,4), 
Ic(5,4)=1. 


YAXIS(5,8), ZAXIS (5,8) 


Figure 3. Input offset points around a ship station 


XAXIS (ISTA) = x locaticz of the station being read in. x=0 is the bow. 


Card 2: FORMAT (2F10.4,75) 
YAXIS(ISTA,I) = 
ZAXIS(ISTA, I) 

Ic(ISTA, 1) 


- location of the Ith point 


; 
z - Location of the Ith point 
ne rth 


~~) 7% 


point is considered a chine roint. 


There should be NT2 aava cards number 2 for each station. 

Note for each paralisl ridbody station, NT2(ISTA) is set equal to 999 
and data card 2 is omitted. This causes the program to copy the results of 
the previous station. 


Sample input data are shown in Appendix II. 


OUTPUT 


A sample output lis<irg is shown in Appendix III. The sample is for the 
input shown in Appendix ==. 

The first set of oucpur is a listing of the input control constants. In 
addition, the computed v2ites for the area of the waterplane (AIL) and the full 
beam (BEAMWL) at the NMIZ svation are listed. These values are used in the 
nondinensionalization ai = later point in the program. It should be noted 
that the table of offsets read in on device number 7 is not reprinted as output. 


This was done in order i> snorten the output. 


es iy Yo 


aE ey ae 
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shiz lergith. Since Bol») Goes net. se Whi Poaato mucher, the value. 
only printed once for each frequency. 

The next set is the three-dinensional source strength (o{x)) distribution 
alors the ship length. The magnitude of o(x) is printed in the column labeled 
MAG(SIGMA). The nagnitude of the nondimensitonal olgma and its real and 


ling Lnary parts are printed in tho subsecioat ihrec columns. ‘he rondinensional 


Sir a is aefined as 
v (x) wg 
TAQ) ee 
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the oressure distribution corer cock stevion ushed Ie: in the input follows 
the source Gistribution listing. if NPNS=-O0, no pressure Gistcibutions are 
princsa. The pressure distribution can only be obtained at the ship stations 
given in the input. The y,z coordinates and angle up from the keel are prin-zad 
in tne first three columis. The twe columns under the heading "2-D POTENTIAI" 
correspond to the two terns of ecuation (7). PRRE squals etY% and PUR! gives 
the values of &y,2;x). The magnitude and phase ansle of the nondimersionai 
pres ire are printed in the last two columns. Thc pressure is computed by 
eguacion (13). The phase angles are all relativ: ‘co = wave rode at the bow. 

If NWAVE=1, the wave amplitude along the shin length is computed and 
printed out. The nondimensional wave amplitude is computed by equaticr (14) 
and includes both the incident pius diffracted wive. As witn the pressure, 
the phase angle is relative to a wave node at the how, 

The sectional exciting force distribution is print out under the reuding 
EXCITING FORCE DISTRIBUTICN. The sectional exciting force (13(x)) is con- 
putei using the expression inside the square brackets of equation (15). Note 
that the e~1¥X 4s not included ‘n the expression for £3(x). In the print out, 
the sectional exciting force is nondimensionalized in the following manner 

aches eB e a 

3 ga B(x) /2 
where B(x)/2 is the LOCAL half beam. The use of the local half bean facili- 
tates conparisons between F3 ara purely two-dimensicnal calculations. 

Finally, the total heave force, pitch moment . .ut midship and the bending 


> 


monent at station number ‘MID are printed out. The results are computed usins 
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equations (16), (17), and (18). The printed results are nondimensionalized in 


the following: manner: 


F* = ¥3 
3 PgaLB 

ye = eee 
pgaL2B 

BM*= BM 
PgaL2B 

where 
L = ship length 
B = full beam at the waterline of the NMID station. 
BM = bending moment at NMID station. 


The phase angles are all relative to a wave node at the bow. 


LIMITATIONS OF THE PROGRAH 


As presently written, the program has several limitations of which the 
user should be aware. The limitations of the subprograms which calculate the 
two-dimensional solution are discussed by Troesch (1976€b). When the hull 
section shape is very thin or has areas of high curvature more input points 
are needed. As the number of input points is increased, the run time of the 
program will be greatly increased. Furthermore, there are eigen frequencies 
at which the solution blows up as described by Troesch (1976b). this is a 
result of the use of an integral-equation technique to solve the two-dimensional 
problem. It only occurs at high frequencies and for normal ship operating 
ranges should be no problem. It should be noted that at springing frequencies 
the effects of the eigen frequencies may become apparent. 

The second major limitation has to do with end effects. As with any 
slender-body theory, the results near the ends are of questionable accuracy. 
The proper means of handling large bulbous bows and transon sterns is nov 
obvious. At the bow, the program can not handle bulb sections which protrude 
underwater forward of the fore-perpendicular. The progran can handle normal 
bulb sections which intersect the free surface. Thus, the user must "fair 


out" the protruding section of the bow. In the stern region, it is importar.+ 


Pee 


that the input data, including the lest station, ft fate in the Jongituctnsa: 
direction. In particular, for cruiser sterms, wore the cftssts of dhe Lary 
station are (0.0, 0.0),it has been observed that the predicted wave amplitude 
sometimes shows a marked change from the previous station. This type of 
result ai the stern section should be viewed with caution. 
The final limitation deals with the forward speed results. The theoreti- 

cal analysis (see equation (3)) was carried out for the two speed ranges T=C 
and «t>>.25. The range between these two extremes has not been formally 
analyzed. The program as presently written has a switch in it at T=.25. i 
For 01.25, the zero speed value of C is used, and for t>.25, C is set equal 
te the targe tvalue. This switch can lead to a diseentinuity in the results 
at tay, 

should be pointed out that the theoretical] analysis assumes high fre- 
Guencics (or short wavelengths). The program appears to give reasonable 
results over the entire freauency range. However, in deriving the expressicn 
for the pressure (equation (13)) certain forward speed terns were neglected 
cecause of the high frequency assumption. Under certain combinations of forward 


speed anc wave frequency these terms may become important. 


-15- 


LIST OF SUBROUTINES 


BESINT - initializes <adles of I Bessel functions for use in determining 


the Green's function for the Helmholtz problen. 
BESx  - computes ths = Bessel function for a given argument and order. 


BIK - is a single nrecision subroutine that evaluates the integral of 


Kg(t) for t ranging from zero to X. 


BKLOD - computes the £4 (X) Bessel function minus its 1/X singularity 


for a given arzunent. 


DBS - solves a system of simultaneous equations using back substitution. 
The L-U decenrosition of the coefficient matrix is computed by 
DLUD. 


DEI, NATSEI, FCNMON, =S=?TA - are subroutines that evaluate the real 
exponential irvesral. They are part of e mathematics package 


from Argonne National Laboratory. 


DETCCH: - calculates <ne comzlex exponential integral. 
DLUD - computes the L-U decomposition of the coefficient matrix by 


Gaussian elinination with partial pivoting. 


DQI4, DQL8, DQL12 - evaluate the integral exp(-X) times some function 
of X for X rarging from zero to infinity by a four, eight and 
twelve poin= Laguerre quadrature formula respectively. These 


subroutines are from the IBM-SSP listing. 


FNTL, FNT2, FNT3 - are all function subroutines that are called by the 


DQL routines. 


G - function useé in computing the integral fo/Vx-E dé 
GRFUN ~ evaluates tne complex Green's function that satisfies the 


Helmholtz ezuation in the fluid domain and the linear free 


surface boundary conditions. 


finds ine netl perareters needed. by tha main propper. These 
ape@lude tye noyvmul, the curvature, ang ths arc dem th for unc 
point in suestior. The paraneters ae; devteruined -.fter 2 


Lircular arc is Tittec through these points. 


inser.s Station: alone shiv length so that a cosine station 


spacing results. 


interpolates values along ship length usirg piecewise cubic 
fit. 


determines the rormal to a line drawn between tuu points. 
coppuces the pressure distribution arcund a ship station. 
evaluates the sine and cosine integrals. 


reads ship offsets from device 7 end comatetsarious se on 


quantities. 


computes tne 3-D source strength by solving the 


intogrel equetior. 


Puen ; eb ci ee se 
a complex function which computes the integral j,9(5)/¥x-* de 


SXY finds a that helps define the source potential in GRFUN. 
x 


; is ~iv 
TRAY - computes tegral fe (xye ax 


TWOATA  - computes 2-D potential,excitiug forces and 


results fron TICDIM. 


computes the solution to the #-l voroblen using rroesch's 


method. 


finds the wave amplitude along the ship length. 
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Program Listing of Diffracted Forces Program 
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